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ABSTRACT 

We examine the stability of a low-mass stellar system surrounding a massive 
central object. Examples of such systems include the centers of galaxies or star 
clusters containing a massive black hole, and the Oort comet cloud. If the self- 
gravity of the stellar system is the dominant non-Keplerian force, such systems 
may be subject to slowly growing (secular) lopsided instabilities. Stability to sec- 
ular modes is largely determined by the dependence of the distribution function 
F on angular momentum J. If dF/dJ < at constant energy, all spherical sys- 
tems are secularly stable. If dF/dJ > 0, as is expected if there is a loss cone at 
low angular momentum, all spherical systems in which F = at J = (an empty 
loss cone) are only neutrally stable, and flattened, non-rotating systems are gen- 
erally unstable. These results suggest that secular instabilities may dominate the 
structure and evolution of the stellar systems in the centers of galaxies. 

Subject headings: Galaxy: center — galaxies: kinematics and dynamics — Oort 
cloud 



1. Introduction 

There is strong evidence that most early-type galaxies contain massive black holes at their 
centers (Kormendy 2004), and tantalizing evidence that smaller black holes may be present 
in some globular clusters (van der Marel 2004). These findings motivate the study of stellar 
systems in which the gravitational force is dominated by a central point mass. We call 
these "near-Keplerian" stellar systems. In this paper we focus on spherical near-Keplerian 
systems since these are a natural first approximation for the centers of elliptical galaxies, 
spiral bulges, and globular clusters. 

Two important features shared by many near-Keplerian systems are: (1) Over a large 
range of radii, the dominant non-Keplerian force arises from the self-gravity of the stellar 



system. (2) Stars on orbits with near-zero angular momentum are destroyed — tidally dis- 
rupted or swallowed whole — by the black hole. In combination with two-body relaxation 
due to stellar encounters, this process leads to an equilibrium stellar phase-space density 
or distribution function (hereafter df) that is an increasing function of angular momentum, 
dF/dJ > 0. If the rate at which relaxation or large-scale torques replenish the lost stars is 
slow enough, there will be almost no stars on low angular-momentum orbits (an "empty loss 
cone"; Lightman & Shapiro 1977). 

A quite different near-Keplerian system is the Oort comet cloud, which surrounds the 
Sun at distances of ~ 10 4 AU. The Oort cloud shares the two characteristic features men- 
tioned above: (1) The mass of the cloud, though quite uncertain, is large enough that over 
much of the cloud self-gravity dominates the other important non-Keplerian forces: the 
quadrupole field from the planets on the inside, and the tidal field from the Galactic disk on 
the outside (see §1.3.2 for quantitative estimates). (2) Comets on low angular-momentum 
radial orbits are lost from the Oort cloud by gravitational interactions with the giant planets. 

The goal of this paper is to investigate whether near-Keplerian stellar systems are subject 
to slow or secular instabilities. These instabilities can arise because orbits in Keplerian 
potentials do not precess, so that even the weak self-gravity of a low-mass stellar system 
can systematically modify the collective orientation of eccentric orbits. The existence and 
nature of secular instabilities is independent of the mass M* of the stellar system relative to 
the black-hole mass M, so long as M* <C M and self-gravity is the dominant non-Keplerian 
force; only the growth rate of any instability depends on M*. 



1.1. Timescales 

We first review the relevant timescales for the evolution of spherical near-Keplerian stellar 
systems. Assume that the system is composed of N = M^/m stars having mass m and 
typical orbital radius r. The dynamical time is tdyn ~ (r 3 /GM) 1 / 2 . If the dominant non- 
Keplerian force is due to the self-gravity of the stellar system, then the orbits precess on the 
secular timescale 

M M 

tsec ~ U yn — ~ t dyn — . (1) 

The two-body relaxation time is (Binney & Tremaine 1987) 

r 3/2 M 3/2 M 2 M 2 

The energies of the stars diffuse on a timescale t re iax- 
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Because the Keplerian ellipse traced out by a star is approximately fixed on a timescale 
t sec , the average torque between two stars is approximately constant over this timescale. 
This slow variation in the mutual torques between stars gives rise to enhanced or resonant 
relaxation of the angular momenta (Rauch & Tremaine 1996); the corresponding timescale 
is 

r 3/2 M l/2 M 
trcs ~ „ U9 ~ tdyn • (3) 

GV 2 m m 

For near-Keplerian systems (M* <C M, N 3> 1) these four timescales are well-separated: 

tdyn tsec tres t rc i ax . (4) 



Most investigations of near-Keplerian systems assume that they are dynamically stable. 
In this case, evolution is driven by relaxation on the timescale t rcs or t re i ax . The basis for 
this assumption is that M* <Mso that the self-gravity of the stellar system is small and 
collective effects should be negligible. However, the small self-gravity from the stellar system 
is responsible for the precession of the orbits and hence there can be instabilities due to 
collective interactions that depend on the distribution of the orientations of the orbits. Such 
instabilities should occur on the secular timescale t sec and hence the growth rate would be 
faster than resonant or two-body relaxation, no matter how small the mass of the stellar 
system may be. 



1.2. Relation between secular, Jeans, and radial-orbit instabilities 

The nature of secular instabilities can be illustrated by a heuristic comparison to two other 
instabilities that affect stellar systems: the Jeans instability that is present in infinite homo- 
geneous systems at wavelengths that exceed the Jeans wavelength Aj, and the radial-orbit 
instability that appears in spherical systems with a preponderance of nearly radial orbits. 

Consider an infinite, homogeneous stellar system. We shall consider perturbations that 
are independent of the spatial coordinates y and z, so the response of the system depends 
only on its equilibrium df F(p), where p is the momentum in the x-direction. We shall 
assume that F is an even function of the momentum. As usual we invoke the Jeans swindle 
(Binney & Tremaine 1987) in which we neglect the contribution of the equilibrium density 
to the gravitational field, so the equilibrium Hamiltonian is H (p,x) = \p 2 /rrii where rrii is 
the inertial mass. We now subject the system to a weak gravitational potential &(x,t) = 
&(t) exp(ikx). The Hamiltonian is modified to H{p,x) = H (p,x) + m g $(x,t), where m g 
is the gravitational mass (usually m, = m g , but we shall use the greater generality). The 
resulting perturbation to the df, f(x,p,t) = f(p,t)exp(ikx), is governed by the linearized 
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collisionless Boltzmann equation, 

9/ xJf „ UJp ^ df dfdH dFd<$> 

where {•, •} is the Poisson bracket. In the present context, this becomes 

f + ^J. ikms if = . ( 6 ) 
at rrii dp 

The perturbed density is p(x,t) = p(t)exp(ikx) where p — fdpf, and Poisson's equation 
reads — /c 2( 3> = AnGp. A neutral mode is present if 

,2 a r [°° dpdF [°°dpdF _ 

k = —AnGm-img / — — — = —QnGmim g I — = kj. (7) 

J-oo P dp Jo P dp 

For the usual situation in which m; = m g , a neutral mode is present if the integral is 
negative, which in turn occurs if the df is a decreasing function of p. More detailed analysis 
shows that the neutral mode divides stable modes with wavelengths A < Aj = 2n/kj from 
unstable modes with A > Aj (the Jeans instability). Thus the Jeans instability is present if 
dF/d\p\ < 0. 

This analysis can be applied by analogy to axisymmetric stellar systems. In this analogy, 
the linear momentum p corresponds to the z-component of angular momentum J z , and the 
position x corresponds to the coordinate conjugate to the angular momentum, which is the 
azimuthal angle of the apocenter, uj. Since x = p/rrii, the quantity corresponding to the 
inertial mass m t is an effective moment of inertia I e s = J z /u; since the precession rate io 
depends on the energy E and angular momentum J z of the orbit, the effective moment of 
inertia is generally also a function of E and J z , which can be positive or negative. 

First consider an axisymmetric stellar system without a central mass point or strong 
central density concentration. In this system, nearly radial orbits cross straight through the 
center of the galaxy, so that each near-radial orbit may be considered to have two apocenters 
separated by Acj) ~ n radians. In terms of the radial and azimuthal frequencies k and Q, 
each apocenter precesses at a rate u — ±(fi — where the two signs apply to prograde 
(J z > 0) and retrograde (J z < 0) orbits. The ratio VL/k is near | for radial orbits whose 
pericenter lies within the central core of the potential, and \ < Q/k < 1 for orbits with 
pericenter outside the core. Thus in general uj > for prograde orbits, so I e & = J z /Cj is 
positive; similar arguments show that I e $ is also positive for retrograde orbits. Since the 
gravitational mass m g > 0, equation (7) suggests that neutral and unstable modes can be 
present if dF/d\J z \ < 0. This is the well-known radial-orbit instability that arises in galaxies 
with radially anisotropic dfs (see Palmer 1994 for a review). 
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Finally we consider a near-Keplerian system. As we suggested at the beginning of the 
paper, in many such systems we expect the equilibrium df to be an increasing function of 
angular momentum, dF/d\J z \ > 0. Equation (7) implies that an instability may be present 
when the product I e sm g dF/d\J z \ is negative. Thus we expect that an instability may be 
present if I e s < 0, which occurs if prograde orbits have to < 0. This is often true in disk 
systems and always true in spherical systems (see §3.1). Thus near-Keplerian stellar systems 
may be susceptible to an instability analogous to the Jeans or radial-orbit instability. 

Although the Jeans instability is always present in a homogeneous stellar system, these 
arguments do not prove that axisymmetric systems with dF/d\J z \ > are always unsta- 
ble. The reason is that in an infinite homogeneous system there can be arbitrarily small 
wavenumbers \k\, while in an axisymmetric system the analog to the wavenumber \k\ is the 
azimuthal wavenumber to, which is restricted to integer values. Modes with m = 1 are the 
most likely to be unstable, but proving stability or instability requires more analysis, a task 
we undertake in §§2-3. 

1.3. Two examples of near-Keplerian systems 

To help relate our theoretical discussion to real systems, we briefly review the relevant 
parameters of two familiar near-Keplerian systems, the Galactic center and the Oort cloud. 
Our main goal is to establish that there is a significant radial range in which (1) the potential 
is near-Keplerian; (2) the dominant non-Keplerian force is self-gravity, and (3) the loss cone 
is empty. 

1.3.1. The Galactic center 

We may describe the mass distribution near the center of the Galaxy by the sum of a point 
mass M = 3 x 10 6 M (the black hole) plus an approximately spherical stellar system with 
density 

p*(r) = Ar~ a . (8) 

Here a — 1.8 and the normalization A can be obtained by setting the enclosed stellar mass 
M*(r) = 4tt J r rfrr 2 p*(r) to 1.4 x 1O 6 M at r = 1 pc (Schodel et al. 2003). With these 
parameters, M*(r) = M at r = 1.9 pc (50 arcsec at 8 kpc) so the stellar system is near- 
Keplerian for r < 1 pc. 

We shall measure radii in milliparsecs or mpc (0.026 arcsec at 8 kpc) since this is not 
far from the accuracy of the best near-infrared positional measurements. For circular orbits, 
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the orbital frequency is given by 

fl" 1 = 0.27rf n / p 2 c y, (9) 

which provides a measure of the dynamical time td yn - The apsidal precession rate for nearly 
circular orbits due to the mass M*(r) is given by 

\oj\~ 1 = 3.8 x 10 3 r°;J c y, (10) 

and this provides a measure of the secular time t sec . The timescale for resonant relaxation is 

M (r) 

tres = tdyn— ~ = 100r^ c J (11) 

assuming a stellar mass m = 1 M . The two-body relaxation time is 

t rclax = 3.5 x 10 8 r r ° n J c y, (12) 

using equation (8-71) of Binney & Tremaine (1987), assuming a Coulomb logarithm log A = 
log(M/m) = 15. 

Solar-type stars are disrupted by the black hole if their pericenter distance is < 9 x 
10 12 cm = 0.003 mpc (eq. [1] of Magorrian & Tremaine 1999). The loss cone is empty for 
solar- type stars for r < 1 pc (based on formulae in §3 of Magorrian & Tremaine 1999). 

The general-relativistic precession rate for nearly circular orbits is given by equation 
(115), 

|o; GR |- 1 = 6.3xl0 2 r^ : y. (13) 
Thus the self-gravity of the stellar system dominates the precession for r > 2 mpc. 

These simple estimates suggest that the three conditions given at the start of this section 
are approximately satisfied in the radius interval 3 mpc < r < 1 pc (0.08" < r < 25"), a 
range of over two orders of magnitude. 



1.3.2. The Oort cloud 

The outer radius of the Oort comet cloud is determined by the ejection of comets by gravi- 
tational perturbations from passing stars, and is probably roughly 5 x 10 4 AU in semi-major 
axis. The inner radius is much less certain, since Oort-cloud comets with semi-major axes 
< 2 x 10 4 AU (the "inner" Oort cloud) are not normally visible from Earth, for reasons 
described below. Models of the formation and evolution of the Oort cloud (Duncan, Quinn 
& Tremaine 1988; Levison, Dones & Duncan 2001) suggest that the inner cloud may extend 
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to semi-major axes of ~ 10 3 AU or even less, with a mass that may exceed the mass of the 
observable outer cloud by a large factor. The total Oort cloud mass of ~ 50M e estimated 
by Weissman (1991) is probably as good as any, but the uncertainty is at least a factor of 
five. The recent discovery of Sedna (Brown, Trujillo & Rabinowitz 2004) suggests that the 
inner Oort cloud may contain more mass and extend to smaller radii than usually believed. 

For our purposes it is sufficient to assume a power-law density distribution of the form 
(8); we shall choose a = 2.5 and normalize so that the mass inside 5 x 10 4 AU is 50M e , 
although steeper exponents are possible (e.g. a = 3.5 is estimated by Duncan, Quinn & 
Tremaine 1988). 

The orbital frequency for circular orbits is given by 

IT 1 = 1.6 x W 5 rf 2 y, (14) 

where r 4 = r/10 4 AU. The apsidal precession rate for nearly circular orbits due to the mass 
of the cloud is given by 

= 9.4 x 10 9 r 4 y. (15) 

The energy and angular-momentum relaxation is dominated by interactions with passing 
stars, and the relaxation time is roughly (Duncan, Quinn & Tremaine 1988) 

t rc i ax = 1 x 10 9 r 4 3/2 y; (16) 

since the dominant interactions are with unbound perturbers there is no analog of resonant 
relaxation. 

Comets are removed from the Oort cloud by gravitational interactions with the giant 
planets, which either eject them from the solar system or perturb them into much more 
tightly bound orbits, if their perihelion distance is < 15 AU. The loss cone due to the giant 
planets is empty if the cometary semi-major axis is < 2 x 10 4 AU. Since the Earth is located 
near the center of this loss cone, Oort-cloud comets are only visible if they arrive from the 
region where the loss cone is full, that is, with original semi-major axes >2x 10 4 AU. Thus 
the properties of the Oort cloud inside this radius are not directly accessible to observation. 

The precession rate of a comet on a nearly circular orbit in the ecliptic due to the 
quadrupole moment from the giant planets is 

IcOpil" 1 = 1.9 x 10 14 r 4 7/2 y. (17) 

Assuming the local density in the Galaxy is pc = 0.10M Q pc~ 3 , the precession rate due to 
the Galactic tide is given by 

|c^ Gal |- 1 = 2.2x 10 9 r 4 " 3/2 y. (18) 
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These estimates suggest that the three conditions given at the start of this section are 
approximately satisfied in the radius interval 200 AU < r < 6 x 10 3 AU, a range of more 
than an order of magnitude. The actual range of validity of our assumptions depends on the 
uncertain normalization and radial dependence of the Oort-cloud mass distribution, 



1.4. Results 

The remaining sections of this paper contain analytical results on the linear stability of 
spherical near-Keplerian stellar systems, which we summarize here. 

1. All such systems in which the df is a decreasing function of angular momentum are 
stable to / = 1 secular modes (§3.4). 

2. All such systems in which the df is an increasing function of angular momentum that 
is non-zero at zero angular momentum are also stable to I — 1 secular modes. 

3. Remarkably, all systems in which the df is an increasing function of angular momentum 
that is zero at zero angular momentum (i.e., the loss cone is empty at all energies) are 
neutrally stable to an I — 1 secular mode (§3.5). The spatial form of this mode 
corresponds to a uniform displacement of the stellar system relative to the central 
mass. 

4. To explore the generality of our results, we have also examined flat, non-rotating, 
near-Keplerian systems, which offer a foil to the spherical systems. We find that flat 
systems do not support neutral modes, and that many flat, non-rotating systems with 
empty loss cones are secularly unstable (§4.4). This result suggests that flattening 
often transforms the neutral mode in spherical systems to an unstable mode. 

5. Numerical normal- mode calculations strongly suggest that secular modes with I > 1 
are generally stable (§4.3). 



2. Linear perturbation theory 

2.1. The equilibrium system 

We consider a spherical stellar system with equilibrium df F(v,x). By "spherical" we mean 
that the df is invariant under spatial rotations of the coordinate frame about the center of 
the system, which implies that F = F(r,v r ,Vt) where r is the radius, v r > is the radial 
speed, and v t > is the tangential speed. 
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We introduce action-angle variables (I, w), which are chosen so that J 2 is the total 
angular momentum, I 3 is the ^-component of the angular momentum, and I± — J r + I2 
where J r is the radial action (Tremaine & Weinberg 1984). The allowed values for the 
actions (the action space) are given by 

< h < h, ~h <h< h- (19) 

The conjugate angle W3 is the longitude of the node, that is, the azimuthal angle at which 
the orbital plane intersects the equatorial plane; there are two nodes separated by n, and 
w 3 is usually chosen to be the one at which the orbit has i > (the ascending node). The 
angle w\ is the orbital phase measured from pericenter, that is, w\ = 2n(t — t p )/T r where T r 
is the radial period and t p is the time of pericenter passage. At pericenter, w 2 is the angle 
measured in the orbital plane from the ascending node. The interpretation of these actions 
and angles in a Keplerian potential is described in §3.1. 

If the df is written as a function of the action-angle variables, then according to Jeans's 
theorem it can only depend on the actions I and on the longitude of the node W3 since these 
are the only integrals of motion in a spherical potential. Rotational invariance requires in 
addition that the df is independent of both W3 and the ^-component of angular momentum 
J 3 . ThusF = F(/ 1 ,J 2 ). 



2.2. The linearized collisionless Boltzmann equation 

We now subject the stellar system to a weak perturbing potential <£>. The resulting pertur- 
bation in the df, /, is determined by the linearized collisionless Boltzmann equation [eq. 
(5), except now the phase space is defined by velocity and position (v, x) rather than mo- 
mentum and position (p,x)]. We restrict ourselves to the case in which the perturbation is 
self-consistent, so that the perturbing potential $ is determined by the perturbed df through 
Poisson's equation, 

0(X , () = _ G /WK^. (20) 

J |x — x'| 

We assume for the moment that the mass M is fixed at the origin; the consequences of lifting 
this assumption are addressed in §3. 

We may assume that the time dependence of / is proportional to exp(At). Equation (5) 
can then be written as a linear eigenvalue equation, 

-{/,ffo}-{F,*[/]} = A/, (21) 

where <&[/] is the linear operator on / defined by equation (20). In action- angle variables 
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this reads 

dw dl dl <9w dw dl dw 



2.3. Expansion in spherical harmonics 

Consider a perturbing potential of the form 

* p (yL) = *i m (r)Y lrn (6,<l>). (23) 
The mass density that produces this potential is 

p p (x) = pim(r)Yi m (9,(f)), 

where 

$!m(r) = _ 2m| dr'{r') 2 Pim{r')-^, P i m {r) = J sin 6d6d<t>Y* m {0, </>)p P (x), (25) 

and r< = min(r, r'), r> = max(r, r'). 

A spherical harmonic can be written in action-angle variables as (Tremaine & Weinberg 
1984) 

i 

Y lm (9, 0) = t m ->yi/ jm (P)e tijx+mW3) (26) 
i=-i 

where 

y y = y y (i7r,0) (27) 

is zero unless I — j is even, r l - m ((3) is a rotation matrix (Edmonds 1960), = cos^ 1 1 3 /I 2 is 
the inclination, and x is the angle measured in the orbital plane from the ascending node. 
Since W2 is a measure of the angle between the ascending node and pericenter, the angle \ 
can be written as u>2 + ip where is a function of I±, I2, and the radial phase w\ only (in §3 
we shall identify ip with the true anomaly in Keplerian potentials). Note that yij and r l j m ((3) 
are both real and y^ is non-zero only if I — j is even. Since the radius r is also a function of 
ii, J2, and wi only, and phase space is periodic in wi, the potential (23) can be written in 
action-angle variables as 

00 1 

$ w= E E^^^^)^^!'^) 6 ^^ 1 ^ 2 ^ 3 ^ ( 28 ) 
ii=— 00 j=—i 

where 

= ^ |* dw, e^-^Ur). (29) 
10 



(24) 



An arbitrary df can be expanded as a Fourier series in the angles, 

/(I, w) = ^ /,(/!, / 2 ,/3)e a w (30) 
1 

Using the expansions (28) and (30), the eigenvalue equation (21) can be solved explicitly, 



i 



i+h-h 



yiiA^Wl^ih.h) (. of of 



*® = — xirmS) + W (31) 

This result shows that the perturbed df arising from a potential perturbation whose angular 
dependence is proportional to Yi m (9, 0) must have the form 

i 

f lm (I, w) = £ ^ ljm 9Um(h, h, w 1 )r l jm (f3)e i ^ +m ^; (32) 
j=-i 



here the factor 



is added to simplify later formulae. 



Vio 



We may now show that the density or potential arising from the response f\ m must 
have angular dependence proportional to Yi m (9,<f>). The response potential having angular 
dependence proportional to Y n k(9,4>) is given by equation (25) as 

4nG f r n 
*nk(r) = ~2n~+l J d^y^'^Vimiv'^)^. (34) 

Since phase-space volume is invariant under canonical transformations, we may replace 
dv'dx' by dl'dw' and substitute equations (26) and (32) for Y nk and fi m . Noting that the 
rotation matrices are real, and using the orthogonality relation 

jT sin (3d(3r l jm ((3)r i ; m ((3) = ^y^, (35) 

we obtain 

® nk (r) = - ^ 1)2 S nl S km \Vu\ J dl'j'^dw', e'^' g ljm {I[, F 2 , w',)^. (36) 

Thus the response potential (or density) is zero unless n = I and k = m. In other words a 
perturbing potential oc Y [m produces a response potential with the same angular dependence, 
a consequence of the rotational invariance of the unperturbed df and potential. 

This result also shows that the eigenfunctions of Poisson's equation and the linearized 
collisionless Boltzmann equation (eqs. 20 and 21) can be chosen to have specific angular 
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dependences: the density and potential can be chosen equal to functions of radius times 
Yi m (6,(f)), and the df can be chosen to be a sum of terms that are functions of I±, I 2 and 
Wi times T l - m {f3) exp[i(jw 2 + mwa)] with \j\ < I. Moreover, the eigenvalues A must be 
independent of the orientation of the coordinate system and hence must be independent of 
m. 



3. Near-Keplerian stellar systems 

We consider motion in an unperturbed Hamiltonian H (v, x) = |t> 2 — GM/r + $*(r), where 
M is the mass of the central object and <l?*(r) is the potential due to a spherical stellar system 
centered on M. The mass of the stellar system interior to r is M*(r) = An J Q r dr r 2 p*(r), 
and we assume that M*(r)/M ~ e <C 1, so that the potential is nearly Keplerian. Since the 
equilibrium system is spherical, GM*{r)lr 2 = d&*/dr. 



3.1. Near-Keplerian motion 

Since the potential is nearly Keplerian, the actions can be written approximately in 
terms of the standard Keplerian orbital elements 

ii ~ (GMa) 1/2 , I 2 ~ (GMa^l-e 2 ) 172 , h = I 2 cosP, (37) 

where a is the semi-major axis, e is the eccentricity, and j3 is the inclination. The conjugate 
angles are 

Wi ~ £, W2 — co, u> 3 ~ fi, (38) 

where £ is the mean anomaly, u is the argument of pericenter, and Q is the longitude of the 
ascending node. If x is the angle measured in the orbital plane from the ascending node 
(eq. 26), then ip = x — co is the true anomaly. 

The Hamiltonian can be written 

^-^ + *,. (39) 
The unperturbed frequencies £1 = w = dH /dl are 

fti ~ (GM/a 3 ) 1/2 , fi 2 ~ 0, fi 3 = 0. (40) 

The frequency Q 2 is equal to the time-averaged value of the precession rate uj, and is 
smaller than Qi by a factor of order e ~ M*(r)/M 1. Although small, f2 2 plays a central 
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role in determining the secular stability of near-Keplerian systems. The simplest way to 
determine fl 2 is through Gauss's method, which states that in a spherical, near-Keplerian 
potential 

. (1 - e 2 ) 1/2 d$* . .... 

^ = n —cosil). (41) 

\liae ar 

Thus 

(1-e 2 ) 1 / 2 /d$, \ G(l-e 2 )V 2 / M.(r)cos^ \ 

where (X) = (2n)~ 1 J Q 27r dwiX represents a time average over the orbit. We can rewrite the 
time average as an average over true anomaly: since dwi = Q\dt, and the angular momentum 
I 2 = r 2 dip/dt, 

( X } = 7T J r d^r 2 X = — — — - / #r 2 X (43) 



2ttI 2 J y 2vra 2 (l-e 2 ) 1 /2 
Using this result in equation (42) we have 



n 2 = -^rr- / d^M,(r)cosij. (44) 
TrMe 



o 



We now investigate the sign of Q 2 - We may write 

dijj M*(r) cos i/j = / d^cos^{M,[r(^)]-M*[r(7r-'0)]}. (45) 



JO 



For < ^ < |7r, r(ir — tp) is greater than r(^). Since M*(r) is an increasing function 
of r, the integrand in equation (45) is negative so f2 2 is also negative. Thus the line of 
apsides precesses in the opposite direction to the orbital motion, so long as the precession is 
dominated by the self-gravity of a spherical system. 



3.2. Secular oscillations 

We now examine secular oscillations in near-Keplerian systems; as described in §1.1, 
these have characteristic growth rates or frequencies of order Q 2 ~ ^ 

We shall work in a frame centered on the mass M. Then in addition to the perturbing 
potential $ described by Poisson's equation (20) there is also an indirect potential arising 
from the acceleration of the frame, which is given by 

*W) = ^7*]^ P( ^ (46) 
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where p — f dv f. However, we now show that the indirect potential can be neglected. 
Since we are interested in secular perturbations, the perturbed density p arises from the 
superposition of the time-averaged density along individual orbits. For any Kepler orbit 

x' \ /cost/A 

^) = (— )•>■ (47) 

where ip is the true anomaly and p is a unit vector pointing toward pericenter. This time 
average is easily seen to vanish, from the last expression in equation (43). Thus the indirect 
potential is zero for secular oscillations, in the approximation at which we are working. Terms 
that are of higher order in the small parameter e are responsible for the weak acceleration 
of the central mass that is required so that the center of mass of M and M* remains fixed. 

Equation (31) provides a formal solution to the linearized collisionless Boltzmann equa- 
tion. To solve Poisson's equation, the perturbed df f\ must be of order the perturbed poten- 
tial Wu 2l:j . However, the unperturbed df F(Ii, I2) is of order e, as is the eigenvalue A and the 
precession rate f^, while the frequency Qi is independent of e. Therefore, self-consistency 
requires that only terms with l\ = are important; in other words, the perturbed df may 
be taken to be independent of w\ and only the orbit-averaged component of the perturbing 
potential W$ 2l is important. The same approximations underlie secular perturbation theory 
in celestial mechanics. 

Using the expressions (28) and (32) for the perturbed potential and df, and restricting 
the eigenvalue equation (21) to terms with l\ — yields 

OF 

-ij^gijm + ij\Vlj I W kimgJ^ = *9ljm, (48) 

where gij m is a function of l\ and I2; from equation (29) 

W? jm = ^f d Wl e ij H lm (r) = (<t>, m (r) cos^>5 (49) 
and $j m is related to gij m by equation (36): 

Sim(r) = "(|^Tj2 E M / MMZtoM,® J dwle-W^. (50) 
These equations can be rewritten as 

OF . - r 

2 j'=-l 
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here we have dropped the subscript m on since the eigenvalue equation is independent 
of m. We have also defined the operator 

T l jf (z) = -j0^\yij\\Vif\ J dl[l' 2 dl' 2 z{l[j' 2 ) J d Wl dw[e^-^^- (52) 

here r and ip are functions of I\, I2, and w\, r' and ip' are functions of I[, I 2 , and w[; and as 
usual r> = max(r, r') and r< = min(r, r'). Note that equation (51) implies that qiq = 0. 

The eigenvalue equation (51) is linear in the eigenfrequency A. This is in contrast to the 
usual equation for normal modes of stellar systems, which is nonlinear in the frequency (e.g., 
Kalnajs 1977; Weinberg 1991). Eigenvalue equations that are linear in the frequency have 
simpler analytic properties and are easier to solve numerically. For disk systems, Polyachenko 
(2004) has stressed that such equations arise whenever the response is dominated by a single 
resonance, that is, by one integer triple 1 in the Fourier expansion in the actions. The result 
(51) is slightly more general, since several resonances with different values of I2 = j (but the 
same values of l± — and Z 3 = m) are involved. 



3.3. The Hermitian form of the eigenvalue equation 

Equations (51) and (52) completely describe the linearized secular oscillations of spherical 
near-Keplerian stellar systems. We now show that this eigenvalue equation can be rewritten 
in a Hermitian form, so long as dF/dI 2 has the same sign everywhere. 

Equation (50) determines the potential $z m (r) in terms of the df specified by gij. The 
true anomaly ip and the radius r are both even functions of w±. Moreover \yij\ is even in j. 
Thus the potential depends on the df only through the combination +gi_j. Therefore it 
is natural to rewrite (51) in terms of the combinations 1 

gt(h, h) = \[giAh, h) + 91,-jih, h)] , 9lj (h, h) = \[giAh, h) - gi,-j(h, h)]- (53) 

Since the operator T 1 --, is an even function of j and j', we find 

dF 1 

-ijn 2 g+ + ij— J2 T IM?) = X 9i P -ij^gij = \gty (54) 
2 j'=-i 

Since the second equation implies that g^ = 0, the sum in the first equation can be restricted 



1 This approach is closely related to Antonov's (1960) trick of analyzing stability of spherical systems by 
forming DFs that are even and odd in the velocities. 
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to f 7^ 0. The two equations can be combined into a single eigenvalue equation, 



A 2 ^ = -j 2 S%g$ + fn 2 §- T l jr {gtr) - SlW- ( 55 ) 

We now restrict ourselves to stellar systems in which dF/dI 2 is non-zero and has the 
same sign everywhere. Since Q 2 < (§3.1) we may define an inner product 

[a ' b] ""5 S pZdf/kf ^ (56) 



3^0 



then S l is a Hermitian operator, since 

r \i /" dl\l 2 dl 2 £l 2 *, 

[ a . s '( b >] = £ y laf7a7T a A 



3,j' = -l " " > 

= [b,S'(a)]*, (57) 



where e = sgn(<9F/<9/ 2 )- Note that the domain of S is restricted to vectors a with a_j = aj. 

The second term in equation 3.3 is directly proportional to the mutual gravitational 
potential energy of the phase-space densities defined by a and b. 

Since S l is Hermitian, its eigenvalues are real and its eigenfunctions can be chosen to be 
orthogonal. Thus A 2 is real, and all modes have one of two forms: either pairs of damped and 
growing modes (A 2 > 0) or pairs of oscillatory modes (A 2 < 0). Some or all of the oscillatory 
modes may be singular or van Kampen modes, as we discuss further in §4.2. 

The stability properties of S l can be investigated analytically for I — 1, as shown in the 
following subsections. A numerical investigation of stability for I > 1 is described in §4.3. 



3.4. Lopsided (I = 1) modes 

The simplest disturbances are those with I — 1; as the arguments in §1.2 suggest, these are 
also the most likely to be unstable. For / = 1 the operator S l can be simplified. 

It is useful to consider the inner product [a, S x (a)]. We have 

/ 3 A 1/2 

yi±1 = T [s^) ' ^o = ' ( 58 ) 
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and a_i = a x . Thus 



dI x I 2 dI 2 £L 2 , l2 



i a - sl < a) i = 2 'mr |a ' 



H — / dlil 2 dl 2 a\ / dl' x l' 2 dl' 2 a! x \ dwidw[ cos ^ cos ^' 



3 ./ " " './ 1 z J './ ' 1 ~~ r r 2 

5a + Sp. (59) 



According to equation (42), 



I 2 f , M*(r) cos^ _ J 2 ^i /" , M*(r) cos 



fi 2 = — =— / d Wl — ^ ^ = / dr w r ; (60) 

27rM e> / r 2 vrMe 7 v r r 2 V ; 

in the last expression we have replaced the integral over the angle Wi from to 2n with twice 
the integral over the radius r from pericenter to apocenter, using the relation dw\ = Q,\dr/v r . 
The radial speed v r is a positive function of I±, I 2 and r defined by v r = [2H (Ii, I 2 ) + 
2GM/r — Jl/r 2 ] 1 / 2 . The first term in equation (59) then becomes 

ttM J r 2 J ev r \dF/dI 2 \ 

Using the relation 

dv r GM ( Jf \ GMecosij 



dr r 2 v r \ GMr J r 2 v r 
this expression can be rewritten as 



(62) 



2 f drM , (r) l ( i h^gi VrW? , (63 ) 



a ttGM 2 J * y J dr J e 2 \dF/dI 2 \ 

where the inner integral is over all values of I\ and I 2 for which v r is real at a given radius. 
We may now integrate by parts; the boundary terms vanish since M* = at r = 0, and since 
we may assume that the perturbation described by a\ vanishes at sufficiently large distances. 
Thus 

where p*(r) is the density of the equilibrium stellar system. 

Similar manipulations on the second term of equation (59) yield 



3GM 2 J y ' ' drdr' _ 

where w(r,r') = (rr') 2 r < /r 2 > = r<. Integrating by parts with respect to r and r', assuming 
that the boundary terms vanish for the reasons given in the preceding paragraph, and using 
the result d 2 w(r,r')/drdr' = 3r 2 5(r' — r), we find 

2 



327T6 

5/3 ~ GM 2 



drr 



dl\l 2 dl 2 Qi 

V r d\ 



(66) 
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First consider systems in which the equilibrium df is a decreasing function of angular 
momentum, dF/dh < 0. Then e = sgn(dF/dh) = — 1 and both S a (eq. [64]) and Sp (eq. 
[66]) are negative. Thus [a, S x (a)] < 0, so the operator S 1 is negative, and all of its eigenvalues 
are less than zero. Thus all spherical near-Keplerian stellar systems with dF/dh < are 
stable to I — 1 secular perturbations. 



Determining stability when dF/dh > requires more work. We use Schwarz's inequal- 



ity, 



with 



Then 



J dhdI 2 \A\ 2 J 
e \dF/dI 2 



dhdh \B\ 2 > 

1/2 



dhdl 2 AB , 
B= {Q lVr dF/dh) 1/2 ■ 



I 



dhdh \B\ 



1 



dhdh QiV; 



dF 
dh 



-I 



dh^Fvrli^o + 



dhhdh ^ 



(67) 
(68) 

'-F; (69) 



in the last expression we have integrated by parts and used the facts that Q± is independent 
of h and 

h 



dh 



r 2 v r 



(70) 



The final integral in this expression is simply related to the stellar density p*{r). To show 
this, we change integration variables from (h, h) to (v r ,Vt) where Vt is the tangential speed. 
We have J 2 = rv t , and at constant J 2 and r the differential energy is dE = Qxdh = v r dv r so 

VrdVr ^ 



dh = 



Thus 



and 



i f dhhdhni r, , „ i 

— / F = / dv r v t dv t F = —p*, 

r 2 J v r J in 



dhdh\B\ 2 < —p*, 
An 



(72) 



(73) 



with equality if and only if the df vanishes at zero angular momentum. Schwarz's inequality 
then becomes 



/ 



dhhdh fii 



vJaA > 



Air 
P* 



dhhdh fii 



-V r di 



(74) 



e 2 dF/dh 

Substituting this inequality into equations (64) and (66) yields 

[a.S^a)] =S a + S p < 0. (75) 

Thus the operator S 1 is non-positive, and all of its eigenvalues are less than or equal to zero. 
We conclude that all spherical near-Keplerian stellar systems with dF/dh > are stable or 
neutrally stable to secular I = 1 perturbations. 
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3.5. The neutral I = 1 mode 



The maximum of [a, S x (a)] is zero, which is achieved when F(/ 1 ,7 2 = 0) = 0, 8F/dI 2 > 0, 
and A — B. This bound is achieved if and only if a is the eigenvector with the largest 
eigenvalue, A 2 = 0. Thus all near-Keplerian stellar systems with a df that is an increasing 
function of angular momentum and vanishes at zero angular momentum support a neutral 
I = 1 mode with 

e OF 

«a K rk (76) 

To convert this to a df, we shall assume that m — 0, that is, we orient the coordinate axes 
so that the I = 1 mode is axisymmetric. The rotation matrices are 

rii,o(/3) = ^sin/3, (77) 

and equation (32) yields 

e dF 

/io oc — — sin/3sincj. (78) 
h oh 

The density corresponding to this df is straightforward to determine. The eccentricity 
vector e, which points toward pericenter and has magnitude equal to the eccentricity, is 
defined by 

v x (r x v) - r. (79) 



GM 

The ^-component of the eccentricity vector is given in terms of the orbital elements by 

e ■ z = e sin /3 sine*; (80) 

and in terms of the velocity by 

e ■ z = [v 2 cos 9 — s r s z v r v z ) — cos 9, (81) 

(jrlVl 

where 9 is the usual polar angle, v is the total speed, v r > and v z > are its radial and 
^-components, and s r (and s z ) are ±1 to account for in-and-out (and up-and-down) motion. 
Thus the df (78) yields the density 

dvf w oc / — — -j^(v 2 cos9 - s r s z v r v z ) - cos9 . (82) 

The equilibrium df F(Ii,I 2 ) depends on the velocity only through v r and v t , where v t is 
the tangential speed. If we define x to be the angle between the tangential component 
of the velocity vector and the unit vector (f) in spherical polar coordinates, then s z v z = 
s r v r cos 9 + v t sin x sin 9 and we have 

f dv r v t dv t dx OF [ r 2 . i 

Pl0(X J T 2 Oh lGM^ Vt C _ Wtsmxsm^) - cose/J . (83) 
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The term proportional to sin x integrates to zero; thus, 



„ f dv r v t dv t dF ( rv? \ 
fto0ccos()j /___^- 1 j. ( 84 ) 

Using equation (71) and the relation rv t = I 2 , this becomes 

cos 9 fdhdh^dF (II \ 

Using equation (62) this simplifies to 

d f OF 
Piooccos^— l dI 1 dI 2 Q, 1 ——v r . (86) 
or J dl 2 

We now integrate by parts with respect to I 2 and assume that F(7 1 ,/ 2 = 0) = since this 
is necessary for the existence of a neutral mode; then using the identity (70) 

pio oc cosu— — - / b occoso 1 — — , (87) 

or r l J v r dr 

where the last expression follows from equation (72). Thus the density in the neutral mode 
is obtained by a uniform displacement of the equilibrium stellar density. 

The neutral mode that we have found is reminiscent of the trivial mode that is present 
in isolated stellar systems when the origin of the coordinate system is displaced from the 
center of the equilibrium system; in spherical systems this displacement mode is an I — 1 
neutral mode with density oc cosOdp^/dr, just as in equation (87). However, the mode 
derived here is quite different, for several reasons. First, it represents a displacement of the 
center of the stellar system from the point mass that dominates the potential rather than 
from an arbitrary coordinate origin (recall that we are working in the frame centered on M) . 
Second, although the spatial structure corresponds to a uniform displacement, the phase- 
space structure does not. Third, the neutral displacement mode is present in all isolated 
equilibrium stellar systems, whereas the mode we have found here is present only in spherical 
systems: for example, disk systems do not support such a mode (§4.4). 



4. Discussion 

4.1. A short derivation of the neutral mode 

There is a simpler proof of the presence of a neutral mode in a wide range of near-Keplerian 
spherical systems. 
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We may orient the coordinate system so that any / = 1 potential perturbation is ax- 
isymmetric, $ p (x) = &i (r)Y 10 (9,<t>). According to equation (28) the secular component of 
this potential can be written as 

$ s (x) = i"V!o(^Wi. h)e ijw *. (88) 
j=±i 

Using equations (29), (58), and (77), this can be simplified to 

/ 3 \ 1/2 

$(x)=f— j ($w(r) costp) sin/^sinu;. (89) 
Thus the secular Hamiltonian is 

/ 3 \ 1/2 

H s (h,I 2 ,u) = H (h,I 2 )+ ( —J ($ 10 (r) cos V) sin /3 since*. (90) 

In a neutral mode, both the Hamiltonian and the df are time- independent. Since the 
secular Hamiltonian is time-independent it is an integral of the motion; the other integrals 
are F and I 3 (since the secular Hamiltonian is independent of w\ and w 3 ). Since the df 
is time- independent, Jeans's theorem implies that it can only depend on the integrals of 
motion, so / = f(H s , I ± , I 3 ). Since the perturbing potential is small we may write 

f^f(H ,I u I 3 ) + -^r(-p\ ' ($ 10 (r)cos^)siii/3sinw. (91) 



dH \Att 

The unperturbed df is F(Ii, I 2 ) and this must equal f{Ho, F, J3). Since Hq depends only on 
Ii and I 2 , the explicit dependence of / on I 3 can be dropped and we must have 

OF = d^dHo =n d£_ 

dh dH dI 2 2 dH ' 1 ; 

Thus the perturbed df is 

fio = f(H s ,I 1 )-f(H ,I 1 )=l—J — — sin/3sinw. (93) 

Now let us assume that the perturbed potential arises from a uniform displacement 
of the unperturbed stellar potential <&*(x) by an amount £ in the z-direction. Then the 
perturbed potential is 

1 /2 

$ p (x) = -^cos^ = f 10 («i) where $ 10 (r) = - ^ (94) 
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Combining this result with equation (42) for Q 2 , equation (93) simplifies to 

GM OF 

f w = -£———es\n(3s\nu. (95) 

J-2 

Following the derivation in equations (79)-(87), the corresponding density is 

p p (x) = J dv/ = -£cos0^. (96) 

where as usual we have assumed that the equilibrium df vanishes at zero angular momentum. 
Comparing this result to equation (94) shows that the response density gives rise to the 
perturbing potential, so we have found a self-consistent normal mode. 



4.2. Oscillatory modes 

In this subsection we discuss the oscillatory secular modes of spherical near-Keplerian stellar 
systems, which have uj 2 = —A 2 > 0. Once again, for simplicity we restrict ourselves to 
lopsided modes (1 = 1). Since the operator S 1 is Hermitian, the set of eigenvalues u 2 (the 
spectrum of — S 1 ) lies on the real line, and since — S 1 is non-negative its spectrum lies in the 
interval [0, 00). 

The properties of the normal modes of stellar systems have been investigated by Mathur 
(1990). In general there is a continuous spectrum, one or more segments of the real line that 
are filled with the eigenvalues of singular normal modes analogous to the van Kampen modes 
of plasma physics. The singular modes that are excited by physical perturbations decay by 
Landau damping. There may also be isolated eigenvalues in the gaps of the real line that do 
not belong to the continuous spectrum, and these correspond to oscillating modes that do 
not decay. 

According to Mathur, the continuous spectrum is the set of all frequencies u 2 spanned by 
the function (1 • Q) 2 over the domain of actions for which F(T) is non-zero. For secular I = 1 
modes in spherical systems, this function simplifies to Q 2 . Thus if Q 2 , ranges continuously 
from 0^ lin to fim^ but not outside this range, isolated oscillatory modes must either have 
< uj 2 < Q 2 miQ or uo 2 > ft* „■ 



There are further constraints on the frequencies of isolated modes. A mode with fre- 
quency uj 2 satisfies [a, S*(a)] = — oo 2 [a, a]. Using equations (56) and (59), 

s H^0k\^-^ (97) 

If uo 2 > ft^ax then the right side must be negative (recall that ^2 < 0). However, equation 
(66) shows that Sp has the same sign as e, which is the sign of dF/dl2- Thus isolated 
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oscillatory modes in systems with dF/dI 2 > must have < u 2 < f^ in , while if dF/dI 2 < 
the modes must have uo 2 > tt 2 



max 



Let us examine further the systems with dF/dI 2 > 0. Since \Q 2 \ — *■ as I 2 — > (eq. 
[42]), isolated oscillatory modes can only exist if the stellar system has a minimum angular 
momentum or maximum eccentricity at every energy, i.e., if the loss cone is empty. It is 
straightforward to show that if the stellar system has density p* oc r~ 7 then the precession 
frequency £l 2 oc a 3 / 2-7 at fixed eccentricity. Thus stellar systems with a minimum semi-major 
axis (set, for example, by collisional destruction of stars), an empty loss cone, and 7 < 1.5 can 
have a gap in the continuous spectrum between f2 min and zero, in which isolated oscillatory 
modes may occur. We have searched numerically for such modes in near-Keplerian stellar 
systems with the df 

f(h I 2 ) oc ( Iitogih/hh) if imin < h< ^ and I 2 > hh\ 
\ otherwise. 

this represents a stellar system in which the loss cone is empty whenever the angular mo- 
mentum is less than a fraction h of the angular momentum of a circular orbit of the same 
energy. If J max J min , the corresponding density is roughly a power law, p oc r~ 7 with 
7 = |(3 — b). The eigenfrequencies were computed using the Kalnajs matrix method (Poly- 
achenko & Shukhman 1981; Weinberg 1991; Saha 1991) using the potential basis functions 
$ p (r) = ru p /(l + r) 3 , where u = (r — l)/(r + 1) and p is an integer (Hernquist & Ostriker 
1992; Saha 1993). Typically 20 basis functions were used. 

For example, consider a df with h = 0.6, 6 = 0, J mm = 1 and J max = (30) 1 / 2 (in units 
where GM = 1). There is a gap in the continuous spectrum given by < uo 2 < f2 min = 14.68, 
in which there is an isolated mode at uo 2 = 8.0 to within 3%. 



4.3. Stability for I > 1 

We have no analytic results on stability for I > 1. However, we have searched numerically 
for unstable modes with / = 2 in stellar systems with the df (98), using Goodman's (1988) 
sufficient criterion for instability. In the current context and notation, this condition can be 
stated as: a spherical near-Keplerian stellar system is secularly unstable if there is a trial 
potential and density $/ m (r) and pi m (r), satisfying Poisson's equation (25), such that 

j=i 

This result does not require any restrictions on the df F(Ii, I 2 ). 
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We have explored a variety of trial functions and several values of the parameters 
-^max/^min! b, and h that define the df. While the search was not exhaustive, we found 
no evidence of instability. 



4.4. Secular modes in flat systems 



Zero-thickness near-Keplerian disks provide an instructive contrast to spherical systems, and 
also provide some indication of how the behavior of spherical systems is likely to be modified 
by flattening. Flat systems are less simple than spherical systems, for two main reasons: 
(1) the apsidal precession rate Q 2 is not necessarily negative, as it is for spherical systems 
(§3.1); thus there is no simple definition of an inner product like equation (56) that makes 
the operator S l Hermitian; (2) the relation between the mass distribution M*(R) and the 
potential (R) is more complicated for flat systems than for spherical ones. 

We assume that motion is restricted to a plane described by polar coordinates (-R, </>). 
We may continue to use the actions I\ and I 2 and the conjugate angles w± and w 2 = uj 
(eqs. 37 and 38), except that now the angular momentum I 2 is a signed variable (I 2 > 
for prograde orbits and I 2 < for retrograde orbits), and the argument of pericenter w 2 is 
measured from the zero of azimuth in the direction of increasing azimuth (instead of from 
the ascending node in the direction of orbital motion). The true anomaly ip is measured 
from pericenter in the direction of increasing azimuth rather than the direction of orbital 
motion; thus the azimuthal angle (f> = w 2 + ip. With these definitions, equation (42) becomes 



where both Vt 2 and I 2 may be positive or negative. 

According to Jeans's theorem, the equilibrium df may be written F(I\, I 2 ); since we are 
interested in analogs to spherical systems we shall assume that the df is an even function 
of the velocity, or an even function of I 2 . Most previous studies of secular stability of near- 
Keplerian disks, which date back to the work of Laplace and Lagrange on the stability of 
the solar system (Murray & Dermott 1999), have focused on systems in which all stars orbit 
in the same direction (Sridhar, Syer & Touma 1999; Tremaine 2001). 

We restrict ourselves to lopsided potential perturbations, of the form 



In the secular approximation, we may average this potential over the fast angle w±, to obtain 




(100) 




(101) 




(102) 
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(compare eq. [89]). Following the arguments in equations (90)-(93), the response df is 

f r = — cos.. (103) 

The response surface density is /i r (x) = f dvf r . To evaluate this we use the relation, 
derivable from equation (79), 

a( Rv I a\ , Rs R v R v 4> ■ , f MA\ 

ecosu = cos0 ) GAI Sm ^ ' 

as usual v R = (2E + 2GM/R - J|/_R 2 ) 1/2 > is the radial speed and s R = ±1 is the sign of 
the radial velocity; and is the azimuthal velocity, which can be positive or negative. The 
contributions of the second term in equation (104) to /i r (x) from s R = ±1 average to zero. 
Thus 

, , n 1 f°° , f°° , 8F cos iP) fRvl \ . . 

We convert this into an integral over action space using I 2 = Rv<p and equation (71), 

*« - W /Iffl^^^^- 1 )' (106) 

Now consider the case in which the potential perturbation (101) corresponds to a uniform 
displacement of the unperturbed stellar potential $*(x) by an amount £ in the x-direction. 
Then 

$ (x) = -£cos0^ = $i(i?) cos0 where $i(i?) = -f (107) 
ai? art 

The corresponding surface-density perturbation is 

// p (x) = -£ cos< ^^T> ( 108 ) 

where pL*(R) is the surface density of the unperturbed disk. Substituting this form for <&i(R) 
into equation (106) and using equation (100) to simplify the result, 

2GM£cos0 [°° [ h dhdF( II \ 

Using equation (62) this simplifies further to 

/i r (x) = -^Rcob^J" ' dmi J* (110) 
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For comparison, the density perturbation (108) that gives rise to the perturbing potential 

is 

Mx) = _, C0S 4 / fcr = -2f™^ / ^F; (111) 

this is not the same as the response density /i r (x) of equation (110) and hence a neutral 
mode with this density distribution does not generally exist in flat near-Keplerian systems. 

The densities that we have derived are useful as trial functions to investigate the secular 
stability of such systems using Goodman's (1988) criterion. For disk systems Goodman's 
criterion states that if there exists a trial surface- density or potential perturbation // p (x) or 
$ p (x) with time dependence exp(At), A > 0, such that the response density £f r (x) satisfies 

J rfx$ p (x)[/x p (x) - ^(x)] > 0, (112) 

then the stellar system is unstable. This result requires only that the df is an even function 
of the velocities. 

In our case we let A — > and use equations (110) and (111); thus 

(113) 

We now integrate the term involving OF/ dl2 by parts, assuming that the loss cone is empty 
(more precisely, we assume that F/I 2 is non-singular as I2 — > 0) so that there is no divergence 
at I2 = 0; the other boundary term vanishes since the upper limit to the range of integration 
is set by vr = 0. Using the relations (62) and (70) we obtain 

//• J(h POO rh J j 

rfx$ p (x)[^(x) - A*r(x)] = 2kGM? J dR-j^- J dm x J ^-|F. (114) 

If this expression is positive then the stellar system is secularly unstable (this result does 
not require either that dF/d\l2\ > or that Q2 < 0). The inner integral is positive, but 
in flat systems d§*/dR is not necessarily positive everywhere (in particular, disks with 
zero surface density near the center, as one expects if the loss cone is empty, have = 
$0 — \aR 2 + 0(i? 4 ) near the center, with a > 0). Nevertheless, many disks have d<&*/dR > 
throughout the radial range containing most of the disk mass, and thus are likely to be 
secularly unstable. 

It has long been known that disk systems with counter-rotating stars are prone to lop- 
sided instabilities, but these analyses have focused on isolated stellar disks or disks embedded 
in a massive halo (Zang & Hohl 1978; Araki 1987; Palmer & Papaloizou 1990; Sellwood & 
Merritt 1994), rather than near-Keplerian disks (Touma 2002). The results from this Section 
imply that lopsided instabilities are also common in near-Keplerian disks with little or no 
net rotation. 
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5. Summary 



Near-Keplerian stellar systems can support secular or slow normal modes, in which the 
eigenfrequency is of order i^(M*/M), where tdyn is the dynamical time, M is the mass of 
the central object, and M* is the mass of the stellar system. We have shown that spherical 
near-Keplerian systems in which the df is a decreasing function of angular momentum, or 
an increasing function of angular momentum that is non-zero at zero angular momentum, 
are stable to / = 1 secular modes. Numerical calculations suggest that these systems are 
also generally stable to modes with I > 1. 

Spherical near-Keplerian systems in which the df is an increasing function of angular 
momentum and zero at zero angular momentum (an empty loss cone) also have no unstable 
modes but, remarkably, are all neutrally stable to an I — 1 mode. The spatial form of 
the neutral mode corresponds to a uniform displacement of the stellar system relative to 
the central mass, although the velocity-space perturbation does not. The analogous zero- 
thickness disks are often unstable to / = 1 modes. 

Other authors have investigated lopsided modes of stellar systems, but in quite different 
contexts from the present paper. Weinberg (1994) has reported nearly neutral / = 1 "slosh- 
ing" modes in linear stability analyses of spherical stellar systems, and long-lived oscillations 
in the centers of iV-body models have been reported by several authors (Miller & Smith 
1992; Spurzem & Aarseth 1996). However, these systems do not contain a central massive 
object so the mechanism responsible for the oscillations is probably quite different. Taga 
& lye (1998) have observed oscillations of a massive central object in iV-body simulations; 
however, they argue that rotation of the stellar system is necessary to drive the oscillation, 
so again the cause of the oscillations is likely to be different. 

Our results suggest that the stellar systems found in the centers of galaxies containing a 
black hole are susceptible to lopsided distortions. Of course, the existence of neutral modes 
does not imply instability, only that the system is close to instability in some sense. Our 
analysis has ignored a number of smaller effects, which may modify the neutral mode to one 
that is weakly stable or unstable. These include: 



Relativistic effects. These induce apsidal precession at a rate 

UJGR ~ a^(l-e^y (U5) 

Since relativistic precession is prograde, while precession due to self-gravity is ret- 
rograde, small relativistic corrections increase — l/f^ and hence enhance instability 
according to Goodman's criterion (99). On the other hand, once the relativistic pre- 
cession rate exceeds the precession rate due to self-gravity it promotes secular stability. 
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Precession due to the giant planets is also prograde, so has similar effects on the Oort 
cloud. 

• Flattening. Flattened, non-rotating systems are likely to be less stable than spherical 
ones, since zero-thickness systems are often unstable (§4.4). 

• Rotation. This probably promotes stability, since rotating near-Keplerian disk systems 
of stars on nearly circular orbits are stable (Tremaine 2001). 

• Non-resonant contributions to the eigenvalue equation. These are smaller by of order 
M*/M and in the usual case when dF/dI\ < will promote stability. 

A natural next step is to examine the secular stability of near-Keplerian stellar systems 
using numerical experiments. This is an arena in which there has been surprisingly little 
activity. Most investigations of near-Keplerian systems focus on their long-term evolution, 
using the Fokker-Planck approximation and assuming spherical symmetry (Bahcall & Wolf 
1976; Cohn & Kulsrud 1978; Freitag & Benz 2001). Recent iV-body simulations of near- 
Keplerian systems do not have sufficient resolution to probe the region in which the loss 
cone is empty (Baumgardt, Makino & Ebisuzaki 2004; Preto, Merritt & Spurzem 2004). 

A related question is whether we can extend the linearized analysis of this paper to 
construct self-consistent models of near-Keplerian systems with significant lopsidedness. The 
general properties of orbits in lopsided near-Keplerian potentials were discussed by Sridhar 
& Touma (1999) but they did not construct self-consistent models. 

If near-Keplerian stellar systems are secularly unstable, then standard models of the dis- 
tribution of stars around black holes and estimates of the rate at which black holes consume 
or disrupt stars may be quite misleading. 
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